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Abstract 

We numerically examine the influence of the viscosity on the relaxation process of localized 
clouds in thermally unstable two-phase media, which are locally heated by cosmic ray and cooled by 
radiation. Pulselike stationary solutions of the media are numerically obtained by a shooting method. 
In one-dimensional direct numerical simulations, localized clouds are formed during the two-phase 
separation and sustained extraordinarily. Such long-lived clouds have been recently observed in 
interstellar media. We demonstrate that the balance of the viscosity with a pressure gradient 
remarkably suppresses the evaporation of the clouds and controls the relaxation process. This 
balance fixes the peak pressure of localized structures and then the structure is attracted and 
trapped to one of the pulselike stationary solutions. While the viscosity has been neglected in most 
of previous studies, our study suggests that the precise treatment of the viscosity is necessary to 
discuss the evaporation of the clouds. 

1 Introduction 

Structures, especially spatially localized ones such as interfaces, shocks, and pulses, play a crucial role 
in understanding elementary processes, dynamical nature and even statistics of widespread nonlinear 
phenomena (e.g., [1, 2, 3, 4, 5,£]). In some cases, these structures are represented by stationary or steady, 
even periodic solutions, e.g. heteroclinic orbit that is a front connecting two states and homoclinic orbit 
corresponding to a pulse. 

In this paper, we try to apply this approach to the formation process of interstellar clouds that are 
cold, dense structures surrounded by warm phase in thermally bistable media. This type of structures 
is often observed in nature: intergalactic clouds and plasmas in tokamaks 8, (see also references 
in the work of Meerson [9j). These cold structures have attracted much attention in recent years in 
astrophysical contexts such as star formation processes (e.g., [TUlfTTlfT^ ). 

Researches on interstellar clouds so far have focused on localized structures but have not system- 
atically dealt with stationary solutions and also their role in the formation process of clouds. This 
is partially because such solutions are not stable and do not survive in direct numerical simulations. 
However, unstable solutions can play a crucial role, e.g., traveling waves and unstable periodic orbits 
in turbulent production of channel flow turbulence [HIEIEII^- We, therefore, try to obtain homoclinic 
orbits that are spatially localized stationary solutions of the governing equations and to describe the 
evolution of neutral and thermally bistable interstellar media by them. 

Our main concern is the effect of viscosity on both formation of and saturation to localized structures, 
although the viscous effect has been neglected or at most implicitly included as numerical viscosity 
because of its expected smallness in interstellar scales. In fact, in the studies so far, such structures are 
generally lost during the separation process [TS1[I1|. However, as shown later, by including the viscous 
effect the evaporation time of the localized clouds is lengthened extraordinarily, where they are trapped 
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close to stationary solutions. Physically speaking, the balance between the viscosity and the pressure 
gradient at a higher order controls the saturation process. This might be related with recently observed 
tiny long-lived clouds [T51 [TBI 112] • 

Here, we will shortly review localized structures in interstellar neutral media. For homogeneous and 
stationary cases, radiative-equilibrium states are attained through the balance of radiative effects at 
each point in space: heating due to cosmic ray and radiative cooling [18j . Field et al. |19j found that 
three of such states exist for a range of pressure and only two of them are linearly stable with respect to 
spatial disturbances. If the system is initially in the unstable state, i.e., thermally unstable one (|20j). 
due to the spatial instabilities the system spontaneously separates into regions of stable states called 
the warm neutral media and the cold neutral media. 

For one-dimensional cases, cold structures are transiently formed during the separation process. Some 
of them merge or evaporate, and finally few clouds survive. They are surrounded by fronts connecting 
two stable radiatively equilibrium states. Zel'dovich and Pikel'ner (hereafter ZP) [3T] studied the single- 
front dynamics in terms of steady solutions with one front in a plane-parallel geometry under an isobaric 
condition. The front is a traveling wave but stationary only when pressure takes a certain value called 
saturation pressure. This nature seems to resemble two-phase systems described by complex Ginzburg- 
Landau equations. However, it should be noted that even in the stationary state, the system is sustained 
by thermal flux compensating excesses of local energy budgets. In this sense, the system is essentially 
thermally nonequilibrium. 

From the view of structures, most of the succeeding studies are based on fronts. Elphick et al. |13] 
extended the ZP solutions to a multifront case with an open boundary condition. In their numerical 
simulation, young cold clouds are merged and a uniform state is finally realized. However, rigid boundary 
conditions [T31 changed the destination of cold clouds. A single system-size cold cloud constituted by 
two fronts was sustained. The conservation of mass in the numerical domain governs the final state. 

For two- or three-dimensional cases, localized structures with symmetries as well as plane- wave fronts 
have been investigated. Graham and Langer obtained spherically symmetric steady solutions un- 
der an isobaric condition. The expanding spherical front, i.e., condensing cold cloud, asymptotically 
coincides with the ZP solution as the radius goes to infinity. Nagashima et al. [T7j also calculated 
the three-dimensional spherically symmetric flow by solving the equations numerically with an open 
boundary condition. They obtained the stationary spherical clouds only when the clouds have critical 
radii depending on pressure, otherwise the curvature effect cannot balance the thermal effects. In or- 
der to compare the spherical clouds to relatively long-lived tiny interstellar clouds observed recently, 
they estimated the evaporation and condensation timescales for contracting (evaporating) and expand- 
ing (condensing) spherically-symmetric fronts of the cold clouds respectively, and concluded that the 
estimated timescales are consistent with the observations. 

This paper is organized in the following way. Section [2] provides governing equations and some 
important physical quantities. Pulselike stationary solutions are numerically obtained in Sec. [31 Section 
[D shows that the solutions are confirmed in direct numerical simulations. Section [5] describes the effects 
of the viscosity and pressure gradient in the dynamics of the localized clouds. We present summary, 
discussion, and conclusion in Sec. [HI 



2 Governing equations 

Consider an optically thin neutral media with external heating and radiative cooling. Self-gravitation, 
magnetic field, and other body forces are neglected. The fluids are described through the following 
basic equations which include the laws of mass, momentum, and energy conservation, together with a 
perfect-gas equation of state: 

| + ^(...)^0. (1) 

d d dP 



2 




Figure 1: (Color online) Radiative-equilibrium curve (thick curve) of the Ginzburg-Landau-type 
heating-cooling function (fTC|) in the TP plane. The heating-cooling function F equals zero on the 
curve, a, 5, and c in Eq. (fT6|) are set to 1.0, 2.0, and 1.0, respectively. The two arrows schematically 
denote that a small perturbation on an unstable state grows, and one of the stable two-phase states 
with lower or higher temperature on the radiative-equilibrium curve is finally realized. 
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where p, v, E, and P are density, velocity, total-energy density, and pressure, respectively. The thermal 
conductivity, the specific-heat ratio, and the viscosity coefficient are denoted by 7, and /i, respectively. 
In general the thermal conductivity is a function of temperature but for simplicity we set it as a constant. 
toh is the molecular mass of hydrogen and ks is the Boltzmann constant. £(T, P) is called a heating- 
cooling function, which describes the total of thermal input and output depending on temperature and 
pressure at each point. Note that fluids can be locally heated and cooled in this system since the fluids 
are sufficiently rarefied. Such fluids are common particularly in astrophysical systems. 
We use the following non-dimensional equations derived from the original Eqs. (Il])-([6]): 



P = pT, (10) 

(11) 
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Here Xj = xj/xq, Uj — uj/uq, t — i/(xo/uo), p = p/po, T = T/T^^, P — P/Pq, and £ — C/Cq, which 
are scaled by the corresponding characteristic values denoted by a subscript zero, respectively. 
The characteristic length scale xo is defined as follows: 



This length is called Field length and represents the front width between two phases [20] . 
We set the characteristic velocity ug to be proportional to a sound speed: 



"0 = \ —Tq. (13) 
V "^H 

There are two non-dimensional constants in the non-dimensional equations: 

7 ksp , . 

EP-^J^- (15) 

Pr is called Prandtl number. Ep represents the degree of thermal contribution to dynamics in Eq. ([9]). 
Larger Ep induces slower dynamics in fluids. 

In this paper, the following Ginzburg-Landau-type equation is adopted for the heating-cooling func- 
tion [13 [ll [11]: 

-pC{T,P) = F{T,P) = a{T-l) - b{T - if - clnP. (16) 

We also tried the more realistic but complicated form of the heating-cooling function ([TT]) used in 
astrophysical contexts [13 [H] : 

pC = -nr + n^A{T), (17) 

r = 2.0 X 10"^*^ ergs s~\ (18) 

= 1.0 X 10^ exp[-1.184 X 10^(7 + 1000)] 

+ 1.4 X 10-2^exp(-92/r) cm^. (19) 

It is found that both the stationary solutions and the results of numerical simulations with the heating- 
cooling function (|17p are qualitatively similar to those with the simplified heating-cooling function (|16p 
(see Figs. 12141 and Figs. 1718]) . Therefore, for simplicity and for elucidation of fundamental mechanism, 
the simplest form of the heating-cooling function defined by Eq. p6p is adopted in this work. 

In Fig. [TJ we show a radiative equilibrium curve on which the heating-cooling function F(T^ P) is 
equal to zero. The S shape of the curve is essential to two-phase separation: supposed to a fixed pressure 
in a range with which the equation has three roots. Because two roots with high and low temperatures 
among them satisfy condition (j20p . they are stable against small temperature (density) perturbations, 

^) <0. (20) 

In contrast, the root with an intermediate temperature satisfies condition (|21|) . and thus a small distur- 
bance on the state grows. This instability is called thermal instability [2U] . 

'OF 



>0. (21) 

dTjp 

Fluid particles around the unstable branch in Fig. [1] finally reach to each of stable branches with high 
or low temperatures via heating or cooling processes. 
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Figure 2: (Color online) Density profiles of pulselike stationary solutions for the three values of pressure 
P = 1.1, 1.01, and 1.001. 

3 Pulselike stationary solutions 

First, we will obtain stationary localized solutions of Eqs. (f7|) — (fTT|) . These are homoclinic orbits. 



3.1 Plane-parallel geometry 

Supposed that a fluid is stationary [d / dr = and = 0), the governing equations are simplified in a 
plane-parallel geometry as follows: 

(22) 



dP 
dx 



= 0. 



(23) 



Because the pressure P becomes uniform via Eq. (|23p . the pressure works as a parameter in the tem- 
perature Eq. ([2H). Here for simplicity and for keeping the temperature positive, we set the parameters 
a, 6, and c of as 1, 2 and 1, respectively. However, note that transforming the variables as x = \/a^i 
T - 1 = , and P = P^/''^ the parameters can be eliminated from Eq. (|22p . 

We impose the following boundary conditions in the semi- infinite space {Q <x < oo): 



dT 
dx 







(x = 0, X — > oo). 



(24) 



Since Eq. (j22p is a second-order ordinary differential equation with respect to x, the two boundary 
conditions in Eq. uniquely determine the solution T{x) if it exists. Simultaneously, the values of 
temperature at the boundaries, i.e., Tq = T{x = 0) and T^o — T{x — > c»), are given. We numerically 
solve Eq. with the boundary conditions by means of a shooting method, and we obtain T{x) 
and also Tq and Too for various values of pressure. In Fig. [21 we show the density profiles 'pix) ~ P/T{x) 
of the solutions for three values of pressure P = 1.1, 1.01, and 1.001. By the refiection symmetry of Eq. 

the solutions are extended to a; < 0. These localized cold solutions are called pulselike stationary 
solutions, hereafter. Figure [3] shows the peak density po — P/Tq and the density at infinity poo = P /Too 
as a function of pressure P. Note that the peak densities are always lower than thejiensity of the stable 
equilibrium state and apart from the radiative-equilibrium curve {F = 0) when P > 1. As P — > 1-K, 
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Figure 3: (Color online) Peak densities po and densities at infinity poo of the pulselike stationary 
solutions. The dashed horizontal lines denote the projections of the solution with P = 1.1 and the ZP 
front solution with the saturation pressure P — 1. 
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Figure 4: Pulselike stationary solutions obtained using the astrophysically realistic but complicated 
form (HZl). 
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Figure 5: (a) Peak densities po and densities at infinity poo of the pulselike stationary solutions of Eqs. 
(|25p — (|27p in two-dimensional and three-dimensional cases, respectively. The two-dimensional and three- 
dimensional densities at infinity almost overlap. The two dashed horizontal lines with P > 1 represent 
the projections of the density profiles of some two-dimensional and three-dimensional solutions into the 
pP plane. The horizontal line with P = 1 represents the one-dimensional ZP front solution, (b) Radii 
of the pulselike stationary solutions as a function of pressure P in one-dimensional, two-dimensional 
axisymmetric, and three-dimensional spherically symmetric systems. 



the peak density asymptotically approaches one of the stable radiative equilibriums, and the pulselike 
solution becomes the ZP solution (front solution). 

We also obtained pulselike stationary solutions using the realistic but complicated form of heating- 
cooling function. The solutions are shown in Fig. HI With the respect that the solutions are also 
homoclinic orbits, they are qualitatively similar to those obtained with simplified heating-cooling func- 
tion (fTB|) as shown in Fig. [51 This supports the use of the simplified heating-cooling function p6p . 



3.2 Axisymmetric and spherically-symmetric geometry 

We also obtain pulselike stationary solutions in multidimensional systems as in the one-dimensional case. 
We suppose that the systems are axisymmetric in two-dimensional cases and spherically symmetric in 
three-dimensional cases. These symmetries make the governing equations simpler as follows: 

nT,P) + ^^ + U = 0, (25) 



dr dr^ 
dP 



•r 



= 0, (26) 



where r and d are the radial coordinate and the space dimension, respectively. The boundary conditions 
are imposed as follows: 

df 

-g~=0 (f=0,f^oo). (27) 

These equations are numerically solved by means of the^ shooting method. As in the case of the one- 
dimensional system, we obtain a localized solution p{r) = P/T{r) with peak density po at the origin and 
density at infinity po^ at r —^ oo. The two densities are dependent only on pressure. Figure [5fa) shows 
Po and Poo as a function of pressure P for the two-dimensional and three-dimensional solutions. As 
pressure increases from the saturation pressure P = I, the peak density leaves away from the radiative- 
equilibrium curve. Note that by the curvature effect represented by the second term in the left-hand side 
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Figure 6: (Color online) (a) Density profiles of tlie warm pulselike stationary solutions of Eqs. and 
for the three values of pressure P = 0.9, 0.99, and 0.999. (b) Peak densities of the solutions in 
one-dimensional, two-dimensional axisymmetric, and three-dimensional spherically symmetric systems 
in the pP plane. The two horizontal lines with P < 1 represent the projections of the density profiles 
of some two-dimensional and three-dimensional solutions into the pP plane. The horizontal line with 
P = 1 represents the one-dimensional ZP front solution. 



of Eq. (j25p . the peak density increases along with the radiative-equilibrium curve as pressure increases, 
unlike the one-dimensional solution in Fig. [3] 

We define the size of the localized structures as a radius Fc, which satisfies T(?c) = {Tq + Too)/2. 
Figure mb) shows the radii as a function of pressure. When P 1+, the radii go to infinity. Since 
the curvature term of Eq. (|25p asymptotically approaches zero in the limit, the localized solutions then 
lead to the ZP front solution in a plane-parallel geometry. The dependence of the cloud size on pressure 
qualitatively agrees with the results of Nagashima et al. [T7] . 

3.3 Bubble solutions 

Warm, rarefied, and localized structures surrounded by cold dense media are also numerically obtained 
in the same way as the cold structures obtained in Sec. 13.21 The density profiles of the solutions for 
the three values of pressure P = 0.9, 0.99, and 0.999 are shown in Fig. Elja). These structures are 
realized only when P < 1 . Figure [D^b) shows the peak densities po and densities at infinity poo as a 
function of pressure P. As in the case of the cold structures, the peak densities are apart from the 
radiative-equilibrium curve when pressure gets smaller than the saturation pressure P = 1. In the 
limit that P ^ 1— , the solutions approach the one-dimensional ZP front solution corresponding to a 
plane-parallel and localized warm structure of an infinity size. These localized warm solutions are called 
" bubbles" [HI ^ and have been frequently referred in astrophysical contexts [5] . 

4 relaxation process of localized cold structures 

We investigate the role of the cold pulselike stationary solutions in the formation of localized structures 
by direct numerical simulations. Equations ([7)l- (|lip are numerically solved in a plane-parallel geometry 
with 8192 grids and a periodic boundary condition. We use a high-accuracy pseudospectral method 
partially dealiased and a fourth-order Runge-Kutta method. As initial conditions, pressure is set to the 
saturation pressure and velocity is zero, that is, P = 1 and u — 0. Density is a small random disturbance 
consisting of only low wave-number modes {k < 16) added to a homogeneous state p — 1 that is on the 
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Figure 7: (Color online) The two-phase separation from the random small density disturbance around 
p — 1- The density profiles at r = 0, 1500 and 165000 are shown. The times correspond to the initial, 
transient, and quasistationary state, respectively. The two horizontal lines represent the densities of the 
states of the radiative equilibrium with the saturation pressure P = 1. 
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Figure 8: (Color online) Formation process of localized structures in direct numerical simulation with 
the astrophysically realistic heating-cooling function (fT7| . The density profiles at the three times t = 0, 
5.0 X 10^ year, and 5.0 x 10'' year are plotted. The two horizontal lines represent the densities of the 
states of the radiative equilibrium at the saturation pressure. 
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Figure 9: (Color online) Density and pressure profiles at r = 165000 in the quasistationary state. The 
pulselike stationary solution with the pressure value at the coordinate where the localized structure has 
a peak density is superposed. 



unstable branch of the radiative-equilibrium curve (see Fig. [Ij. We set non-dimensional constants Ep 
and Pr to 1.0 x 10^ and 1.0, respectively. Koyama and Inutsuka0 used Ep ~ 64 in a realistic situation 
so that Ep — 100 seems to be a feasible value. The specific-heat ratio 7 is set to 5/3. 

Figure [7] shows a formation process of localized structures. The initial density disturbance rapidly 
grows and evolves into warm and cold regions due to a thermal instability. Some of the cold regions 
merge successively and some of them evaporate in the evolution. Several localized cold structures finally 
survive in warm media, and then these structures persist for a long time. Because most of the kinetic 
energy in the whole system is lost during the two-phase separation, the cold structures are almost 
stationary. 

For comparison, we also have numerically calculated the two-phase separation with the realistic but 
complicated form of heating-cooling function (fT7|) and original astrophysical Eqs. ([I])—©. Some of 
the results are shown in Fig. [H The relaxation process is qualitatively similar to that with simplified 
heating-cooling function (|16p shown in Fig. [71 Some of the cold regions merge and evaporate, and some 
localized structures finally survive in the process. 

To examine the decaying process of localized structures in detail, we focus on two typical localized 
cold structures. One with a longer lifetime is seen around x ~ 60 in Fig. [71 The peak density of this 
localized structure does not yet attain the value of the density in the radiative equilibrium {F = 0). 
However, as seen in Fig. [SI the localized cold structure is quite similar in form to a pulselike stationary 
solution with the value of pressure at the peak of the localized structure. The values of pressure at 
the peak of localized structures are called peak pressure. Note that in Fig. [51 the value of pressure is 
subtracted by 1.000 042 5 so that the pressure variation around the peak of the localized structure is 
negligible and does not affect the estimation of the peak pressure. 

The peak density of the localized structures decreases slowly and monotonically by evaporation and 
they finally disappear. We estimate the characteristic time of the decaying process, i.e., the evaporation 
time of the localized structures. In Fig. [TUl we show the evolution of the peak density of the localized 
structures obtained in two runs under the same conditions except for the initial condition: one with a 
longer lifetime is called LSI and the other LS2 hereafter. The evolution is separated into two periods: 
the fast transient one and the slow decay one. We call the latter a quasistationary state. The transient 
period finishes at around r = 2.0 x 10'^. This time corresponds to a sound crossing time defined as a 
time during which a sound wave passes through the system. 

The decay curves of the peak density after attaining the quasistationary state are well fitted by a 

^H. Koyama and S. Inutsuka, e-print ^arXiv: astro-ph /0605528 , 
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Figure 10: (Color online) Peak density evolutions of the localized structure at x ~ 60 in Fig. [7] and 
the other localized structure derived by another run. Inset: the peak density evolutions magnified in 
the quasistationary state 4.5 x 10^ < r < 1.1 x 10^. The curve with the larger peak density is fitted by 
the double-exponent function defined as ^ exp{— i? exp[(r — To)/rc]} with A — 3.27, B — 1.24 x 10^^, 
and the relaxation time Tc — 2.91 x 10^. The curve with the smaller peak density is fitted by the above 
double-exponent function with A ^ 2.95, B ^ 1.07 x 10"^, and = 2.21 x 10"*. 



double-exponential function Aexp{— i?cxp[(T — to)/tc]} as shown in the inset of Fig. [TOl the values of 
the fitting parameters of LSI are A = 3.27, B = 1.24 x 10"^ tq = 45000, and n = 2.91 x 10^ for 
4.5 X IQ-* < T < 1.6 X 10^ those of LS2 for 4.5 x 10" < r < 1.1 x 10^^ are A = 2.95, B = 1.07 x 10"^ 
tq = 45000, and Tc = 2.21 x 10"*. The latter case shows that the late decaying process is faster than 
a double-exponential function, but the characteristic time estimated is shorter than the period of the 
saturation process 1.2 x 10^. In this sense, the characteristic time gives a good estimate of the evaporation 
time of the localized structure. 

It should be noted that the relaxation times Tc are much longer than the sound crossing time. In 
Sec. m we discuss the details of this decaying process of the peak density. 



5 persistence of pulselike clouds with viscosity 

Why can the localized cold structures be sustained extremely longer than expected? Here we focus 
on viscosity, although the viscosity has been omitted in most of the previous studies. Since we have 
performed numerical simulations with a high-accuracy spectral scheme, we can deal with the viscosity 
more precisely than the previous studies. 

Figure [TT] shows the viscosity term and the pressure-gradient term of the kinetic equation ([S]) around 
the localized structure in the quasistationary state r = 165000. While the pressure is almost uniform 
because of the mixing by sound waves, a weak pressure gradient is maintained around localized structures 
(see also Fig.[S]). Moreover, the weak pressure gradient almost coincides with the viscosity. This balance 
between the viscosity and the pressure gradient causes the persistence of the localized structure. To see 
this in detail, we rewrite the basic equations ([7]) and ((51) in a plane-parallel geometry as follows: 

^ = (28) 

, 3^2 f4:ld 4 1 dp\ /' 7 - 1 Pr 5? 



Dt ' 4^^ \3pdx 3pdx)[ 7 Ep dx ' ' ^^^^ 
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Figure 11: (Color online) The viscosity term and the pressure-gradient term of the kinetic Eq. 
the quasistationary state with Prandtl number Pr = 1.0 for LSI at r = 165000. 
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Figure 12: (Color online) The time evolution of the viscous stress a. The curve is fitted by an exponential 
function: ctq exp((r — ro)/Tco-), where tq = 70000, ctq = 6.20 x 10~*, and Tea — 3.63 x 10^. 
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Figure 13: (Color online) Comparison of localized structures and corresponding pulselike stationary 
solutions in the relaxation process shown in Fig.[Tni (a) LSI and (b) LS2. The value of pressure of the 
pulselike solutions is set to the peak pressure of the localized structure at each time. 



where D/Dt = d/dr + vd/dx is a Lagrangian derivative. 5 ~ 4:/3{dv/dx) is a viscous stress, which 
represents the degree of compressibility. 

Figure [T^] shows that the viscous stress a exponentially grows in the quasistationary state. Hence, 
the Lagrangian time derivative Da/Dr can be written as aa, where the growth rate is a '--^ 10~^. The 
very small growth rate is just caused by the small difference between the viscosity and the pressure 
gradient in Eq. ([29]) when the nonlinear term in the left-hand side is negligible. In fact a ^ aa ^ 
since the viscous stress a is the order of 10~^ as seen in Fig. [T^l The slow exponential growth of the 
viscous stress finally leads to the double-exponential decay with a long relaxation time via Eq. (pS)) . 

In this relaxation process, as shown in Fig. [T31 the localized structures remain close to pulselike 
stationary solutions. The balance fixes the peak pressure of the localized structure. Then the structure 
is attracted and trapped to one of the pulselike stationary solutions with the corresponding pressure 
value that seems to behave similar to a saddlelike fixed point in the phase space. The deviation of the 
localized structure from the pulselike stationary solution induces both the small pressure gradient and 
viscous effect (see Fig. [TT|) . Finally the two effects induced balance to each other and then this balance 
suppresses the induction of flow that transfers heat and material from localized structures. 

The one-dimensional Eqs. (^5]) and ([^ can be extended to multidimensional equations as follows: 



Dp 

Dt 



-(divu)p. 



(30) 



_D(div u) duj dui 
Dt dxi dxj 

f i A _ ( tzIEl^ _ ^ (31) 

\ p dxi (P' dxi / I 7 Ep dxj dxi j 

These equations suggest that the balance between the viscosity and the pressure gradient plays a crucial 
role in the relaxation of the multidimensional localized structures through the divergence of velocity. 
Therefore, the viscosity should be considered as well as the pressure gradient regardless of the dimension 
of the system. 
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Figure 14: (Color online) Evolutions of peak pressure and peak density of the localized structures LSI 
and LS2 shown in Fig. (TUl The equilibrium curve and the pulselike stationary solution illustrated in 
Fig. [3] are also drawn. The points on the evolution curves are plotted every 5000 from r — 5000 to 
T — 165000 for LSI and from t — 5000 to t = 135000 for LS2. Inset: the evolutions of peak pressure 
and peak density magnified around the times in which the localized structures are trapped near the 
corresponding pulselike stationary solutions. 

6 Discussion and concluding remarks 

We have investigated localized cold structures in thermally unstable two-phase fluids. The two-phase sep- 
aration is induced by thermal instability via local heating and cooling described by the Ginzburg-Landau- 
type heating-cooling function. We have numerically obtained pulselike stationary solutions of which the 
peak densities are not in the states of the radiative equilibrium in one-dimensional, two-dimensional 
axisymmetric, and three-dimensional spherically symmetric systems using a shooting method. Our so- 
lutions differ from the solutions obtained in the previous studies which consider bistable media and fronts 
connecting the two stable states [131 HH HT]- The size of the solutions depends only on pressure, and 
the relationship between the size and the pressure qualitatively coincides with the results of Nagashima 
et al. PT] , 

We have carried out direct numerical simulations of the thermally unstable fluids with a small random 
density perturbation using a high-accuracy spectral method. During the two-phase separation, many 
cold clouds are formed, and some of them merge or evaporate. Finally, several cold clouds are left in 
the warm media and persist for a long time. The pulselike stationary solutions obtained by the shooting 
method are actually observed in the direct numerical simulations. Such small and long-lived clouds have 
been observed in interstellar media and might be related with the long-lived localized clouds we have 
found [HI [13 [HI- 

We have shown that the viscosity balances with the pressure gradient in the quasistationary state 
at higher order. The balance remarkably suppresses the evaporation of the localized clouds. Note that 
the viscous and the pressure-gradient effects are very weak but play a crucial role in the persistence of 
the structures. In contrast, most of the previous studies have neglected the viscous effect (e.g., Refs. 
[T31[T11[33[T71[35]). Our results strongly suggest that the viscosity changes the relaxation process and 
suppresses the evaporation of the cold clouds. We, therefore, expect that multidimensional clouds, which 
are also pulselike as the one-dimensional clouds, have a long relaxation time because of the viscosity. 

Although the pulselike stationary solutions are obtained under constant pressure, the pressure is a 
dependent variable in direct numerical simulations. However, the balance between the viscosity and the 
pressure gradient tends to keep the pressure nearly constant (pressure saturation). In this situation, 
one of the pulselike stationary solutions is selected and attracts the localized structure. This solution 
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behaves similar to a saddlelike fixed point in phase space. Figure [M] shows the evolution of peak values 
of the localized structures LSI and LS2 in density-pressure space. The peak pressure and density of 
the localized structure approach those of the pulselike stationary solution. After a long stay around 
there, it quickly leaves from there to the lower equilibrium state. Note that the leaving process of 
LSI is not seen in Fig. dH In this final escaping process, the profile of the localized structure quickly 
deviates from the pulselike stationary solution with the same pressure value as the peak pressure of the 
localized structure (see Fig. I13|) . By this picture of the relaxation process, the relaxation time might be 
dependent on the approaching process along with a stable manifold of the pulselike stationary solution. 
When similar to the case of LS2, the pressure saturation is not sufficient, that is, localized structures 
are not close enough to the stable manifold of one of the pulselike stationary solutions, the localized 
structures seem to approach a part of a manifold formed by the pulselike stationary solutions rather 
than a single one, i.e., a fixed point. Anyway, the dynamical system approach will help us understand 
the formation process of localized clouds. 

In the future, we plan to perform direct numerical simulations in two-dimensional axisymmetric and 
three-dimensional spherically symmetric systems in order to confirm the persistence of the multidimen- 
sional localized structures with the viscosity. In addition, we will employ two-dimensional full numerical 
simulations in order to investigate multidimensional effects such as interface dynamics. 
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